function [W_end,S,t] = cruise(ac,W_start,ISA)
global oswald

% STEP BY STEP PROCEDURE FOR DET. AIRCRAFT RANGE
% REF.: ROSKAM A & P - P.527
% ACCURATE DETERMINATION OF RANGE BY NUMERICAL INTEGRATION OF SPECIFIC
% RANGE AND ENDURANCE
% SI USAGE

%% INPUT DATA


%%% Internal function variables %%%%

%INITIAL CRUISE WEIGHT DEFINITION
Wcrz_beg=W_start/9.81; %Aircraft cruise initial weight - [kg]
W_empty=ac.weight.empty/9.81; %Aircraft empty weight - [kg]
[Cd0]=calc_Cd0_aircraft(ac,ac.cruise.altitude,ac.cruise.speed); %Aircraft parasite drag coefficient
crz_range_tgt=ac.cruise.range/1000; % Aircraft cruise range [km]

% Internal variables
S = ac.geo.wing.S;
a = ac.geo.wing.ar;
e = oswald;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% REPLACE
g = 9.81;
h = ac.cruise.altitude; %Altitude in [m]


V_tas=(50:5:200); %TAS array - used in ROSKAM methodology [m/s]

%% Load Engine data for Mach & Altitude, full throttle assumed
for k=1:1:size(V_tas,2)
throttle=ac.eng.th_cruise; %Cruise throttle assumed
speed=V_tas(k);
% % % % % % [engine] = engine_model(ac); %Access engine function
[TSFC,T,Wf_dot] = engine_model(ac,h,speed,throttle,ISA); %Access engine function
TSFC_crz(k)=(TSFC); %Engine thrust specific fuel consumption - [kg/N.s]
end

%% REQUIRED THRUST


i=1;
for W=Wcrz_beg:-((Wcrz_beg-W_empty)/25):W_empty %Weight array

%REQUIRED POWER

temp=0;  %Delta T from ISA
flag=1 ; %-> temp means teh absolute temperature in K
%Atmospheric properties, SI units
[p,rho,mu,T,sound,delta_isa]=isa_cos(h,temp,flag);

q=0.5*rho.*V_tas.^2; %Dinamic pressure array
CL=(W*g)./(S*q); %Lift coeficient array
CD=Cd0+((CL.^2)/(pi*a*e)); %Drag coeficient array

% REQUIRED POWER - REQUIRED THRUST

T_req=(CD.*q.*S); %Required Thrust [N];
P_req=(CD.*q.*S).*V_tas; %Required power [W];
SHP_req=P_req/735; %Required power [SHP];

%% CRUISE FUEL FLOW

    for k=1:1:size(V_tas,2)
    Wf_dot_cruise(k)=TSFC_crz(k)*T_req(k); %Total engines fuel flow - [kg/s]
    end

%% SPECIFIC RANGE CALCULATION

SR(i,:)=((V_tas)./(Wf_dot_cruise))/1000; %Specific range = [km/kg]

V_tas_p(i,:)=V_tas; %V_TAS plot array

%% SR x V graph, CONSTANT ALTITUDE - RANGE DEFINITION

W_store(i)=W; %Weight plot array, store current weight value

i=i+1; % Next weight

end

%figure(1)
%for j=1:size(SR,1)
%plot(V_tas_p(j,:),SR(j,:))
%hold on
%end
%xlabel('TAS [m/s]');ylabel('SR [km/kg]');
%hold off
%grid;

%% SR x Weight integration
error=1e6; %Initial cruise range error (%)
factor=0.99; %Weight end cruise = factor * Weight cruise start - initial value

while abs(error)>0.01 %0.01% of range deviation

Wcrz_end=Wcrz_beg*factor; %End of cruise weight estimative
kk=1;

    for W=Wcrz_end:50:Wcrz_beg %For cruise weight defined array
    SR_p(kk)=interp2(V_tas,W_store,SR,ac.cruise.speed,W); %Interp SR
    W_p(kk)=W; %Store current plot weight
    kk=kk+1;
    end

crz_range=trapz(W_p,SR_p); %Cruise range numerical integration - [km]
error=100*((crz_range-crz_range_tgt)/crz_range_tgt); %Error percentage
    if error>0
        factor=factor+(0.0001*abs(error)); %Correct factor
    end
    if error<0
      factor=factor-(0.0001*abs(error)); %Correct factor
    end
% error %%%%%%Display error
end

%figure(2)
%plot(V_tas,SR)
%xlabel('TAS [m/s]');ylabel('SR [km/kg]');
%legend(num2str(W_p'))
%grid;
%figure(3)
%plot(W_p,SR_p)
%xlabel('W [kg]');ylabel('SR [km/kg]');
% grid;

%% Function Output

W_end=Wcrz_end*9.81; %End of cruise weight - [N]
S=crz_range*1000; %Cruise range - [m]
t=S/ac.cruise.speed; %Cruise time - [s]

%% Extra Outputs
t_hour=t/3600; %Cruise time - [h]
Wf_consumed=(Wcrz_beg-Wcrz_end); %Total amount of fuel consumed - [kg]
SR_mean=mean(SR_p); %Specific range - [km/kg]
SR_mean_imperial=SR_mean/4.083; %Specific range - [lb/nm]
